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SUMMARY 


The  work  reported  herein  was  initiated  in  response  to  the  interest  and 
need  expressed  by  the  Aerotheraodynamics  personnel  at  the  Naval  Weapons  Center, 
China  Lake,  California.  This  report  contains  the  results  obtained  in  the  first 
phase  of  the  project  on  aero-heating  calculations  where  only  the  hemisphere- 
cylinder  configuration  was  studied.  Other  configurations  will  be  studied  in 
the  future. 


This  effort  was  supported  by  Naval  Air  Systems  Command  (NAVAIR)  and 
executed  for  the  Naval  Weapons  Center  under  the  Strike  Warfare  Weaponry 
Technology  Block  Program  under  Work  Request  N00019-78-WR-81079,  Airtask 
A03W-03P2/0Q8B/7F32-3000-0Q0  (appropriation  178  1319.1981).  This  airtask 
'provides  continued  exploratory  development  in  the  air  superiority  and  air-to- 
surface  mission  areas.  Mr.  W.  C.  Volz,  AIR  320C,  was  the  cognizant  NAVAIR 
Technology  Administrator. 


Useful  discussions  with  C.  F.  Markarian,  W.  R.  Compton  and  B.  Ryan  in  the 
course  of  this  work  are  gratefully  acknowledged. 
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1.  INTRODUCTION 


Excessive  thermal  stress  is  known  to  be  the  main  cause  of  mechanical  failure 
of  most  long-range,  high-speed  missiles.  Therefore,  in  the  design  of  aerodynamic 
configurations  for  advanced  tactical  missiles,  the  ability  to  accurately  predict 
the  aerodynamic  heating  is  essential. 

The  objective  of  this  project  is  to  develop  efficient  engineering  procedures 
for  predicting  the  aerodynamic  heating  on  missile  dome  configurations  related  to 
the  Navy’s  tactical  missiles.  The  heating  under  consideration  consists  of 
convection  from  the  boundary  layer  on  the  surface  of  the  dome  and  conduction 
inside  the  structure;  the  two  modes  of  heat  transfer  are  coupled  at  the  surface 
of  the  dome.  In  certain  high  speed  missiles,  plastic  radomes  are  used  which  may 
be  subject  to  ablation/charring  during  their  long  range  flight.  Therefore,  the 
capability  of  treating  conduction  problems  with  phase  transitions  is  also 
desirable. 

In  the  present  development  of  such  engineering  methods  of  aerodynamic  heating 
calculations,  the  basic  guidelines  are  simplicity  and  accuracy.  The  simplicity 
is  dictated,  to  a large  extent,  by  economic  considerations  because  the  method  is 
to  be  used  extensively  for  parametric  studies.  The  starting  point  of  the  present 
work  is  the  procedure  of  aerodynamic  heating  calculations  currently  used  at  the 
Naval  Weapons  Center-*-  (NWC).  The  central  part  of  this  computation  procedure  is 
the  finite  difference  code  for  two-dimensional  conduction  calculations,  known  as 
the  SINDA  code.  The  necessary  boundary  conditions  for  the  conduction  problem  are 
provided  by  the  empirical  solution  of  the  convective  heating  associated  oith  the 
boundary  layer  flow  over  the  dome  configurations.  However,  the  NWC  procedure 
appears  inadequate  in  certain  applications,  and  the  need  for  its  Improvement  has 
become  obvious. 

In  compliance  with  the  aforementioned  guidelines,  the  efforts  in  the  present 
task  are  primarily  directed  towards  developing,  or  improving,  approximate  (or 
empirical)  techniques  for  aerodynamic  heating  calculations.  We  pursue  the 
objective  in  two  directions,  one  being  to  use  the  existing  technology  in  the 
area  and  the  other  being  to  develop  new  technology  to  meet  future  needs. 

In  this  report,  some  preliminary  results  of  the  effort  are  reported. 

Section  2 contains  a description  of  a proposed  scheme  for  laminar  convective 
heating  calculations  and  some  comparisons  with  the  scheme  currently  used  by  NWC. 
While  the  new  scheme  is  based  on  existing  formulas  and  is  still  empirical  in 
nature,  the  accuracy  seems  significantly  better  than  that  of  the  current  NWC 
scheme.  An  evaluation  of  the  available  methods  for  transitional  and  turbulent 
heating  calculations  is  also  included  in  this  section.  In  section  3,  a simple 
integral  method  for  transient  heat  conduction  calculations  is  presented  along 
with  its  applications  to  a class  of  idealized  ablation  problems.  This  method 
is  presented  here  as  a potential  candidate  for  use  in  the  heat  conduction 


1.  Compton,  W.  R. , "Aerodynamic  Heating  of  Spherically-Tipped  Cylinders,  Cones 
and  Ogives,  Using  the  General  Thermal  Analyzer  SINDA,"  NWC  'IN  4061-172, 

Jun  1974. 
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calculations  in  the  future.  The  results  of  this  section  represent  the  efforts 
in  developing  new  technology  for  aerodynamic  heating  calculations.  Finally, 
some  concluding  remarks  and  future  plans  are  stated  in  section  4. 

The  three  principal  dome  configurations  used  in  tactical  missiles  are  shown 
in  Fig.  1.  In  this  paper,  results  of  convective  heat  transfer  are  presented 
only  for  the  hemisphere-cylinder  configuration.  Application  of  the  procedure  to 
other  configurations  is  not  believed  to  cause  basic  difficulties  as  long  as  the 
surface  pressure  formulas  can  be  satisfactorily  verified,  and  results  of  the 
application  will  be  reported  in  the  future. 

2 . AERODYNAMIC  HEATING  PREDICTION  METHODS 


2 . 1 Surface  Pressure  Distribution 


The  starting  point  for  predicting  aerodynamic  heating  on  a missile  dome  is 
the  prediction  of  the  surface  pressure  distribution  on  the  dome,  which  is 
subsequently  used  to  compute  the  local  flow  properties  at  the  edge  of  the 
boundary  layer.  The  dome  geometries  considered  are  all  blunt  bodies  with  detached 
shock  waves  at  supersonic  speeds,  and  local  flow  properties  may  be  computed  from 
a specified  surface  pressure  distribution  since  the  total  pressure  behind  the 
normal  shock  wave  at  the  stagnation  point  is  related  directly  to  the  flight  Mach 
number.  In  Ref.  2 an  empirical  equation  developed  by  Andrews^  is  shown  to  give 
accurate  predictions  of  surface  pressure  on  hemispheres  and  other  bodies  of 
revolution  at  Mach  numbers  above  .75.  This  empirical  equation  is  a modification 
of  Newtonian  flow  theory,  and  has  the  following  form: 


C 

P 


cos  0 + 


h 

Si 


.78 

M2 3-27 


sin0cos0  - 


.95  sin0 


exp [2. 235 (M  -1)] 


(2.1) 


where 


C 

P 


C 

po 


M_ 


P 


P 


CO 


2(P  “ Pj 

local  pressure  coefficient  = — = 

Y*T  p 

stagnation  point  pressure  coefficient 
total  length  of  nosetip,  ft 
freestream  Mach  number 

2 

local  surface  pressure,  lbf/ft 

2 

freestream  static  pressure,  lbf/ft 


R^  = body  radius  at  end  of  nosetip,  ft 
Y = ratio  of  specific  heats  for  air  (y  = 1.4) 

0 = local  angle  between  dome  surface  normal  and  axis  of  symmetry. 


2.  Isaacson,  L.  K.  and  Jones,  J.  W. , "Prediction  Techniques  for  Pressure  and 
Heat  Transfer  Distributions  over  Bodies  of  Revolution  in  High  Subsonic  to 
Low  Supersonic  Flight,"  NWC  TP-4570,  Nov  1968. 

3.  Andrews,  J.  S.,  "Steady  State  Airload  Distribution  on  a Hammerhead  Shaped 
Payload  of  a Multistage  Vehicle  at  Transonic  Speeds,"  Boeing  Co.  Rpt. 
D2-22947-1,  Feb  1964. 
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The  accuracy  of  equation  (2.1)  has  been  further  tested  in  the  present  investiga- 
tion by  comparison  with  several  additional  sources  of  pressure  data  for  the 
hemisphere-cylinder  configuration4>5»6,7,8  including  recent  theoretical  and 
experimental  results  from  AEDC  (Ref.  4).  This  comparison  is  shown  in  Figures  2a 
and  2b  and  it  can  be  seen  that  the  empirical  equation  compares  favorably  with  the 
other  sources  of  data.  A simple  exponential  relationship  was  shown  in  Ref.  2 to 
give  good  results  downstream  of  the  shoulder.  A comparison  of  predictions  from 
this  type  of  empirical  relation  with  several  sources  of  data  is  shown  in  Figure  3. 
The  empirical  equation  used  is  as  follows: 


°P  ' (Cp>90"  exp 


" wn,  1 

" 2(m-.4)J 


(2.2) 


where 


(V90° 


pressure  coefficient  at  the  shoulder  (from  Eq.  (1)) 
m = 1.2,  for  M < 1.2 

CO 

= M , for  M>  1.2 


Ax  = axial  distance,  measured  from  shoulder 
Note  that  on  the  cylinder  approaches  zero  as  Mm  increases. 

The  comparisons  made  indicate  that  equations  (2.1)  and  (2.2)  can  be  used  to 
obtain  reasonably  accurate  predictions  of  surface  pressure  on  hemisphere-cylinders 
over  the  Mach  number  range  of  interest  (1  to  5)  in  missile  dome  aerodynamic 
heating  calculations.  In  using  equation  (2.1),  however,  it  was  found  that  a 
minor  modification  was  necessary  to  obtain  satisfactory  results  in  the  region 
near  the  stagnation  point. 

The  local  velocity  at  the  edge  of  the  boundary  layer  computed  using  equation 
(2.1)  and  one-dimensional  flow  relations^  was  found  to  be  linear  with  distance 


4. 


5. 


6. 


7. 


8. 


9. 


Hsieh,  T.,  "Flow-Field  Study  About  a Hemisphere-Cylinder  in  the  Transonic  and 
Low  Supersonic  Mach  Number  Range,"  AEDC  TR-75-114,  Nov  1975. 

Baer,  A.  L.,  "Pressure  Distributions  on  a Hemisphere  Cylinder  at  Supersonic 
and  Hypersonic  Mach  Numbers,"  AEDE  TN-61-96>  Aug  1961.  . 

Katz,  J.  R. , "Pressure  and  Wave  Drag  Coefficients  for  Hemispheres,  Hemisphere- 
Cones  and  Hemisphere-Ogives,"  NAVORD  Rpt.  5849,  Mar  1958. 

Stine,  H.  A.  and  Wanloss,  K. , "Theoretical  and  Experimental  Investigation  of 
Aerodynamic-Heating  and  Isothermal  Heat-Transfer  Parameters  on  a Hemispherical 
Nose  with  Laminar  Boundary  Layer  at  Supersonic  Mach  Numbers,"  NACA-TN-3344, 

Dec  1954. 

Morrison,  A.  M. , et.al.,  "Handbook  of  Inviscid  Sphere-Cone  Flow  Fields  and 
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along  the  surface  (or  angle  0)  over  most  of  the  hemisphere,  as  shown  in  Figure  4a. 
This  is  in  agreement  with  results  from  other  investigations^.  Close  to  the 
stagnation  point  however,  the  velocity  is  nonlinear  and  at  Mach  numbers  above 
about  1.5  a value  of  zero  velocity  occurs  at  some  distance  from  the  stagnation 
point.  This  behavior  is  due  to  the  fact  that  the  empirical  equation  gives 
erroneous  values  of  Cp  close  to  the  stagnation  point,  and  values  of  Cp  greater 
than  C occur  at  Mach  numbers  above  1.5,  as  shown  in  Figure  5.  It  should  be 
pmax 

noted  that  according  to  one-dimensional  flow  theory  (Ref.  9): 


where 


dU 
e 

de 


1 dp 

p U d9 
e e 


(2.3) 


Ue  = local  velocity  at  edge  of  boundary  layer,  ft /sec 
Pg  = local  density  at  edge  of  boundary  layer,  lbf-sec^/ft^ 

Since  Ue  -*■  0 as  0 -*■  0,  dUe/d0  can  remain  finite  only  if  dp/d0  0 (and 

hence,  dCp/d0  0).  This  condition  is  satisfied  by  the  pure  Newtonian  equation, 

but  not  by  the  empirical  terms  which  have  been  added  in  equation  (2.1).  In  the 
present  investigation  this  difficulty  was  circumvented  by  projecting  the  lines 
of  Ue/U„  versus  0 to  0 = 0 and  shifting  each  line  by  a constant  value  of  AUg/U^ 

so  that  Ug/U^  * 0 at  0 = 0.  The  correction  factors  used  are  shown  in  Figure  4b. 

Local  pressures  and  pressure  coefficients  consistent  with  the  adjusted  values  of 
local  velocity  were  computed  and  used  in  the  aerodynamic  heating  calculations. 

The  adjusted  values  of  Cp  have  been  indicated  in  Figures  2 and  3 by  dashed  lines, 
and  it  can  be  seen  that  the  changes  in  local  Cp  due  to  these  adjustments  were 
small. 


The  calculation  of  aerodynamic  heating  at  the  stagnation  point  requires  that 
the  velocity  gradient  at  the  stagnation  point  be  known.  Values  of  the  non- 
dimensional  velocity  gradient  0D/Um*  were  evaluated  from  the  curves  of  Ve/Vm 
versus  0 shown  in  Figure  3 and  have  been  compared  with  results  from  another 
investigation^  in  Figure  6.  The  agreement  is  quite  good. 

2.2  Aerodynamic  Heating  Calculations 

Aerodynamic  heating  rates  calculated  by  the  methods  currently  in  use  at  the 
Naval  Weapons  Center  (NWC)l  have  been  compared  with  predictions  from  several 
other  methods  by  the  use  of  a specific  example.  The  example  chosen  was  that  of 
a 3-inch  diameter  hemisphere-cylinder  flying  at  Mach  numbers  ranging  from  1 to  5 
at  an  altitude  of  approximately  20,000  ft.  The  conditions  for  the  calculations 
are  listed  in  Figure  7, 


10. 


* 0 


Korobkin,  I.,  "Laminar  Heat  Transfer  Characteristics  of  a Hemisphere  for  the 
Mach  Number  Range  1.9  to  4.9,"  NAVORD  Rpt.  3841,  Oct  1954. 
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ds7o  ■ R (df)o  ■ I Wvo’  1/Sec 


* distance  measured  along  surface,  ft 
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2.2.1  Laminar  Heating 


The  current  prediction,  method1  of  NWC  is  first  briefly  described  in  the 
following.  The  method  utilizes  the  stagnation  point  heating  equation  of  Fay  and 
Riddell 


r'o  - ■U': 


'ooV4  vt° 


-vVK 


(2. A) 


where 


2 

Q = stagnation  no:.nt  heating  rate.  Btu/ft  sec 
o 

Pr  = Frandtl  number  for  air 

3 

,=  air  density  evaluated  aL  stagnation  pressure  and  temperature,  lbm/ft 

Uo  = viscosity  of  air  evaluated  at  su-gnation  temperature,  lbm/ft-sec 

cr  air  density  evaluated  at  sti'gue*ion  pressure  and  wall  temperature, 

’A'  lb /ft 3 


p,  " vj.scosi  y of  air  evaluated  at  wall  temperature,  lbm/ft-sec 

0,  = heat  capacity  of  air  at  constant  pressure  = .24  Btu/lbm  °R,  assuming 
ideal  gas  properties 

Tq  = stagnation  temperature,  °R 

Tw  * wall  temperature,  °R 


/due\ 

\T^~J0  ~ stagnation  point  velocity  gradient,  1/sec 

Heating  rates  downstream  of  the  stagnation  point  on  the  hemispherical  dome  are 
computed  from  the  equation  developed  by  Lees^  based  on  a modified  Newtonian 
pressure  distribution: 


wh  -.re 


(2.5) 


11. 


12. 


Q » lv-;ai  heating  rate,  Btu/ft  -sec 
f (0)  = 


l - e2  - 


9sin48  1 - cos49 

oT  *■ 


8 


[»2- 

. „ - 


6sin29  + 


1 - cos29 


] 


Fay,  J.  A.  and  Riddell,  R.  F.,  "Theory  of  Stagnation  Point  Heat  Transfer  to 
Dissociated  Air,"  Journal  of  Aerospace  Sciences.  Vol.  25.  No.  2.  do.  71-85. 
Feb  1958.  

Lees,  L.,  "Laminar  Heat  Transfer  Over  Blunt-Nose  Bodies  at  Hypersonic 
Flight  Speeds,"  Jet  Propulsion.  Vol.  26,  No.  4,  pp.  259-269,  Apr  1956. 
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On  the  cylinder,  the  current  method  employs  a serai-empirical  relationship 
developed  by  Kays^  for  aerodynamic  heating  on  axisymmetric,  bodies.  The  follow- 
ing equation  is  based  on  the  tabulated  values  given  in  Ref.  13  with  the  Prandtl 
number  dependence  incorporated  for  easy  usage  in  later  discussions 


_ /T  \-.12  /T  \-.< 
c.  . 33  lawl  f w i 

* ‘ Pr2«  /R - lT3  / \ aw/ 


(2.6a) 


where 


p U C (T  - T ) 
e e p aw  w 


/S  r2(p  U )1,87  dS 
o e e 

2 . rt  \ ^ 

p r (p  U ) 
e e e 


(2.6b) 


Tflw  * the  adiabatic  wall  temperature,  °R 
Tg  = the  local  temperature  at  the  edge  of  the  boundary  layer,  °R 
P£  = the  viscosity  evaluated  at  the  local  edge  temperature,  lbm/ft-sec 
r * the  local  body  radius,  ft 

Heating  rates  computed  from  equations  (2.5)  and  (2.6a)  are  matched  at  the 
hemisphere-cylinder  junction  and  a constant  pressure  and  local  edge  velocity  are 
assumed  on  the  cylinder. 


Although  equation  (2.6a)  is  only  utilized  on  the  cylinder  in  the  method 
currently  used,  its  range  of  application  Includes  axisymmetric  bodies  of  general 
shape  and  according  to  Kays^3  it  may  also  be  used  at  the  stagnation  point.  At 
the  stagnation  point  equation  (2.6a)  reduces  to  the  following: 


Q * .728  Pr"2/3  (T  /T  )"*08  (p  p )'5  C (T 
o w o'  0 0 p o 

then,  comparing  equations  (2.4)  and  (2.6c): 


- vV(ar)„ 


(2.6c) 


(VFay-Riddell  ^ ,76  Pr 


/VwY1 

\poV 


13.  Kays,  W.  M. , Convective  Heat  and  Mass  Transfer,  McGraw  Hill,  1966 
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t 


t 


and, 


^VlCays 


A comparison  of  Q0  values  from  equations  (2.4)  and  (2.6 r)  over  a range  of 
possible  flight  conditions  is  given  in  the  table  below: 

Pr  = .7,  T = -12°F 

00 

(VKays  ~ ^o^Fay-Riddell 


M 

00 

T /T 
w 0 

^o^Fay-Riddell 

1 

.9 

-1.2% 

2 

.6 

- .4% 

5 

.35 

+2.5% 

In  this  comparison,  the  Sutherland  viscosity  law  for  air  has  been  used  to 
evaluate  Pw/v0.  It  can  be  seen  that  equations  (2.4)  and  (2.6c)  are  in  good 
agreement  over  this  range  of  conditions. 

Predictions  of  laminar  aerodynamic  heating  by  the  method  currently  used 
have  been  compared  with  predictions  obtained  from  a finite  difference  computer 
code1^  and  with  predictions  based  on  the  consistent  use  of  equation  (2.6a)  over 
the  entire  hemisphere-cylinder,  utilizing  local  flow  properties  computed  by  the 
procedure  outlined  irv  section  2.1.  The  comparison  of  laminar  heating  predictions 
by  these  three  methods  is  shown  in  Figures  8a  and  8b.  The  results  based  on  the 
consistent  use  of  equation  (2.6a)  have  been  labeled  "proposed  procedure"  in  the 
figure.  Comparisons  are  shown  for  Mach  number  2 and  above,  since  aerodynamic 
heating  on  missile  domes  is  considered  relatively  unimportant  below  Mach  2.  A 
slight  inconsistency  in  this  comparison  should  be  mentioned.  In  the  finite 
difference  calculations  a value  of  .72  was  used  for  Pr,  whereas  a value  of  .7 
was  used  In  the  current  method  and  proposed  method  calculations.  This 
inconsistency  introduced  roughly  a 2 percent  difference  in  the  results,  which  is 
hard  to  detect  on  the  scale  to  which  the  results  have  been  plotted  in  Figure  8. 

The  finite  difference  method  is  considered  to  be  accurate  for  laminar  heating 
and  is  used  here  as  a standard  of  comparison.  It  is  seen  that  the  predictions 
based  on  the  presently  proposed  scheme  are  in  good  agreement  with  the  finite 
difference  method,  particularly  at  the  lower  Mach  numbers,  and  offer  a significant 
Improvement  over  the  method  currently  used.  Recent  flight  test  data  suggest  that 
boundary  layer  transition  may  not  occur  on  most  IR  domes  at  high  altitudes,  and 
consequently,  the  prediction  of  laminar  heating  rates  should  be  of  practical 
importance  and  interest. 


14.  Cebeci,  T.,  Smith,  A.  M.  0.  and  Wang,  L.  C. , "A  Finite-Difference  Method  for 
Calculating  Compressible  Laminar  and  Turbulent  Boundary  Layers,"  McDonnell 
Douglas  Aircraft  Co.  Rpt.  DAC-67131,  Mar  1969. 
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2.2.2  Turbulent  Heating 


For  turbulent  heating,  the  method  currently  used  employs  a semi-empirical 
relationship  developed  by  Kaysl3: 


where 


St  « 


♦ 0295 

Pr,4(R  Y1 
eT 


(p  U )r1,25dS 
e e 


Ue  r_ 


1.25 


(2.7a) 


(2.7b) 


Predictions  by  the  current  method  have  been  compared  with  predictions 
obtained  from  the  finite  difference  computer  code-^  and  with  predictions 
obtained  from  an  equation  developed  by  Vaglio-Laurin^  which  is  similar  to  that 
of  Kays: 


St* 


(2.8a) 


where 


St' 


p U C (T  - 
e e p o 


V 


/S  (p  U )ii  r1*25  dS 
o e e e 

2 1.25 

r 


(2.8b) 


The  comparison  of  turbulent  heating  predictions  for  Mach  numbers  2 and  5 is 
shown  in  Figure  9.  In  this  case  the  finite  difference  method  cannot  be  used  as 
a standard  since  empirical  relations  concerning  turbulent  transport  must  be 
utilized.  The  three  methods  compared  appear  to  be  in  reasonably  good  agreement, 
especially  near  the  shoulder  of  the  hemisphere  and  on  the  cylinder.  Further 
assessment  of  the  turbulent  prediction  methods  must  involve  comparisons  with 
experimental  data. 

2.2.3  Transitional  Heating 


Downstream  of  the  point  where  boundary  layer  transition  begins  skin  friction 
and  heat  transfer  rates  increase  rapidly  from  their  laminar  values  to  values 


15.  Vaglio-Laurin,  R.,  "Turbulent  Heat  Transfer  on  Blunt-Nosed  Bodies  in  Two- 
Dimensional  and  General  Three-Dimensional  Hypersonic  Flow,"  Journal  of 
Aerospace  Sciences,  Vol.  27,  No.  1,  Jan  1960. 
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anproaching  those  for  a completely  turbulent  boundary  layer.  The  rapid  increase 
of  heating  rates  in  the  region  of  boundary  layer  transition  is  an  important 
factor  in  determining  thermal  stresses  in  missile  domes,  making  it  desirable  to 
predict  heating  in  this  region  as  accurately  as  possible.  Several  procedures 
exist  for  predicting  transitional  skin  friction  and/or  heating  rates,  all  of 
which  are  more  or  less  empirical  and  are  generally  based  on  experimental  data 
from  flat  plate  experiments.  In  the  following  discussion,  it  is  assumed  that 
the  location  where  transition  starts,  S*,  is  given. 

a.  Persh/NWC  Method 

Persh^  developed  a method  for  computing  boundary  layer  growth  in  the 
transition  region  based  upon  the  following  equation  for  the  skin  friction: 


0 


f 


constant 


(2.9) 


A 


< 


The  growth  of  the  boundary  layer  in  the  transition  region  is  computed  using 
equation  (2.9)  and  the  boundary  layer  momentum  equation: 

d9  _£f  e L+2>  1 dUe  ■ 1 dPe  , 1 dr]  (2 

dS  ~ 2 6 (H+2)  U dT+p  dS  rds| 

Lee  j 

where 

= local  skin  friction  coefficient 

0 = boundary  layer  momentum  thickness 

Re  = Reynolds  number  based  on  momentum  thickness 
0 

H * boundary  layer  shape  parameter. 

Subscript  "T"  denotes  the  fully  turbulent  value 

The  constant  in  equation  (2.9)  is  determined  by  the  condition  C.  * Cf  at  the 

L 

start  of  transition,  where  Cf  is  the  laminar  value  of  the  skin  friction 

tL 

coefficient.  Accordingly,  equation  (2.9)  may  be  rewritten  as  follows: 


16.  Persh,  J.,  "A  Procedure  for  Calculating  the  Boundary^Layer  Development  in 
the  Region  of  Transition  from  Laminar  to  Turbulent  Flow,"  NAVORD  Rpt.  4438, 
Mar  1957. 
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or, 


where 


Superscript  * denotes  values  at  the  start  of  transition 


(2.11a) 


(2.11b) 


Persh  compared  predictions  from  this  procedure  with  experimental  data  for 
incompressible  and  compressible  flow  over  flat  plates  and  obtained  reasonably 
good  agreement  in  both  cases.  In  order  to  apply  the  method  to  the  more  general 
case  of  boundary  layer  flow  on  axisymmetric  bodies  some  method  of  evaluating  the 
boundary- layer  shape  parameter,  H is  required.  Persh  developed  a method  for 
predicting  H based  on  a correlation  of  velocity  profile  data  for  transitional 
boundary  layers  on  flat  plates.  The  Involved  nature  of  the  required  calculations 
makes  the  Persh  method  rather  inconvenient  to  use  except  when  combined  with 
procedures  which  compute  laminar  and  turbulent  boundary  layer  growth  by  numerical 
integration  of  the  integral  boundary  layer  equations. 


In  using  the  Persh  method  for  predicting  transitional  heating  it  is  assumed 
that  the  parameter  F has  the  same  value  for  heating  as  for  skin  friction: 


F » 


Cf  ' Cf 
T 

C„  - C„ 


Tj 


StT  - St 
StT  ' StL 


The  current  method  used  at  NWC'L  to  predict  transitional  heating  rates  is  a 
modification  of  the  Persh  method,  where  the  parameter  F is  related  to  the 
equivalent  turbulent  flat  plate  length  Reynolds  number  as  follows: 


F = 


(2.11c) 


The  integration  used  to  obtain  Re^,  (equation  2.7b)  is  modified  as  follows: 
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so  that  at  S*: 


I*  + / (p  U Jr1* 25  dS 


Rem  = 


S* 


e e 


*e  r 


1.25 


(2.7c) 


* 

ReT  = 


I* 


p*  r* 
e 


1.25 


(2. 7d) 


The  quantity  I*  is  determined  by  requiring  Reg  to  be  the  same  for  the 
turbulent  case  as  for  the  laminar  case  at  S*.  In  the  NWC  procedure: 


for  the  laminar  boundary  layer 
for  the  turbulent  boundary  layer 
combining  the  above  gives: 

R*t  - 36.9 


Reg  = .664  /Re^ 

(2.12a) 

Re0  « .037  (ReT)4/5 

(2.12b) 

)5/8 

(2.12c) 

and  I*  follows  from  eq.  (2.7d). 

At  S*,  the  value  of  Re-p  from  equation  (2.12c)  is  used  in  equation  (2.7a)  to 
compute  St£.  At  locations  downstream  of  S*  values  of  Stx  are  computed  from 
equations  (2.7a)  and  (2.7c).  Thus,  the  turbulent  values  of  St  used  in  the 
definition  of  F are  based  on  an  effective  starting  point  for  the  turbulent 
boundary  layer  which  is  downstream  of  S = 0.  In  the  Persh  method  an  identical 
procedure  is  used  to  determine  the  turbulent  values  of  the  skin  friction 
coefficient. 

b . Spot  Theory 

Experimental  investigations  of  transition  on  flat  plates  indicate  that  the 
transition  region  is  characterized  by  the  intermittent  appearance  of  turbulent 
spots,  which  grow  as  they  move  downstream  and  finally  merge  to  form  the  turbulent 
boundary  layer.  Emmons-*-'  developed  a theory  for  predicting  skin  friction  in  the 
transition  region,  based  on  the  growth  of  turbulent  spots,  involving  an  inter- 
mittency  factor  defined  as  the  fraction  of  time  a given  location  is  occupied  by 
turbulent  spots.  All  averaged  properties  such  as  skin  friction  and  heating 
rate  adjust  smoothly  from  laminar  to  turbulent  values  as  the  intermittency  factor 
increases  from  0 to  1.  Dhawan  and  Narasimha-*-®  using  flat  plate  data  found  that 


17.  Emmons,  H.  W. , "The  Laminar-Turbulent  Transition  in  a Boundary  Layer  - 

Part  I,"  Journal  of  Aerospace  Sciences,  Vol.  18,  No.  7,  pp.  490-498,  Jul  195L 

18.  Dhawan,  S.  and  Narasimha,  R.,  "Some  Properties  of  Boundary  Layer  Flow  During 
the  Transition  from  Laminar  to  Turbulent  Motion,"  Journal  of  Fluid  Mechanics, 
Vol.  3,  1957-58. 
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the  origin  of  turbulent  spots  takes  place  very  nearly  along  a single  line  across 
the  flow  located  at  the  beginning  of  transition.  They  found  a universal  inter- 
mittency  distribution  for  flat  plate  flows  and  successfully  predicted  the  skin 
friction  and  velocity  profiles  in  the  transition  region.  Chen  and  Tyson-^ 
extended  Emmons  spot  theory  to  include  the  case  of  boundary  layers^  on  blunt 
bodies.  The  parameter  F is  related  to  the  intermittency  factor,  y>  from  the 
turbulent  spot  theory  as  follows: 


F = 1 - y 


Chen  and Thyson  give  the  following  relationships  for  1 - y, 


or 


for  the  flat  plate: 


(2.13) 


(2.14a) 


where  G = the  spot  formation  rate  parameter,  l/ft-sec 
R = hemisphere  radius,  ft 


(2.14b) 


Chen  and  Thyson  give  the  following  relation  for  the  dimensionless  spot  formation 
rate  parameter: 


where 


(Res)2 

r*2*68a2 

ee 


* 

Re, 


* * * 

p u s 

e e 

y* 


A - 60+4.68  M 


*1.92 


(2.15a) 


(2.15b) 


Equations  (2.15a)  and  (2.15b)  were  obtained  from  a correlation  of  the  extent  of 
transition  Reynolds  number  |e^ g as  a function  of  the  transition  Reynolds  number 
Reg  and  edge  Mach  number,  Me,  using  data  from  flat  plate  experiments. 


19.'  Chen,  K.  K.  and  Thyson,  N.  A.,  "Extension  of  Emmons’  Spot  Theory  to  Flows 
on  Blunt  Bodies,"  AIAA  Journal,  Vol.  9,  No.  5,  pp. 821-825,  May  1971. 
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c.  Comparison 

A comparison  was  made  between  values  of  the  parameter  F computed  from  the 
Persh  method,  the  current  NWC  method  and  the  Chen  and  Thyson  method  for  the  case 
of  incompressible  flow  on  a flat  plate — this  being  the  only  case  where  results 
can  be  obtained  in  a straightforward  manner  from  the  Persh  method.  The  defini- 
tion of  F in  terms  of  the  skin  friction  coefficient  was  used  in  the  Persh  and 
NWC  methods  and  the  following  relations  for  incompressible  flat  plate  boundary 
layer  flow  were  utilized: 

. . . , . Cf  .322  .220  x 

laminar  boundary  layer:  «—  = = — — (2.16a) 

\fRe~  0 

\ x 


turbulent  boundary  layer:  ~ = _ _•  (2.16b) 

(Rex)1/5  Re01/4 

where  Rex  = the  length  Reynolds  number  for  flat  plate  flow 

The  results  of  comparisons  made  for  values  of  the  momentum  thickness 
Reynolds  number  at  transition  of  200  and  400  are  shown  in  Figures  10a  and  10b. 

In  the  figures,  the  parameter  F has  been  plotted  against  the  flat  plate  length 

Reynolds  number  measured  from  the  beginning  of  transition,  ARe  . 

X 

Of  the  three  methods  compared,  the  Chen  and  Thyson  method  predicts  the 
shortest  extent  of  transition  while  the  NWC  method  predicts  the  longest  extent. 

The  Persh  method  and  Chen  and  Thyson  method  are  in  fairly  good  agreement  over 
the  initial  part  of  the  transition  region  at  the  transition  Reynolds  number  of 
400.  Plotting  ARex  on  a log  scale  as  in  Figures  10a  and  10b  makes  it  difficult 
to  compare  the  gradients  of  skin  friction  predicted  by  the  three  methods.  The 
NWC  method  predicts  much  higher  gradients  of  skin  friction  (and/or  heating)  at 
the  start  of  transition  than  the  Chen  and  Thyson  method,  but  lower  gradients 
toward  the  end  of  transition.  In  the  NWC  method  the  steepest  gradients  occur 
near  the  start  of  transition  whereas  in  the  Chen  and  Thyson  method  they  occur 
toward  the  middle  of  the  transition  region. 

A final  comparison  of  predictions  for  transitional  heating  is  shown  in 
Figure  10c.  Here,  predicted  values  of  the  parameter  F are  shown  for  the 
transitional  boundary  layer  flow  on  a hemisphere  at  Mach  5.  The  conditions 
used  in  the  boundary  layer  calculations  are  those  shown  in  Figure  7 for  a 
Mach  number  of  5,  and  transition  was  assumed  to  start  50  degrees  from  the 
stagnation  point.  Predictions  by  the  Persh  method  have  not  been  included  due  to 
the  complexity  of  the  calculations  required  for  flow  on  a hemisphere.  In 
Figure  10c  predicted  values  of  the  parameter  F from  the  NWC  method  and  Chen  and 
Thyson  method  have  been  plotted  against  the  equivalent  turbulent  flat  plate 
Reynolds  number  measured  from  the  start  of  transition,  ARe^,  obtained  from 
equation  (2.7c). 

Predictions  from  the  Chen  and  Thyson  method  using  both  the  equation  for  a 
flat  plate  (2.14a)  and  the  equation  for  a hemisphere  (2.14b)  are  shown  in 
Figure  10c.  The  difference  between  the  two  curves  is  very  small  indicating 
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there  is  little  gained  by  using  the  more  complicated  function  for  the  hemisphere. 
The  comparison  between  the  NWC  and  Chen  and  Thyson  methods  is  essentially  the 
same  as  the  comparison  for  skin  friction  in  incompressible  flat  plate  flow. 

In  consideration  of  the  comparisons  made  here  it  is  felt  that  the  method 
currently  used  at  NWC  probably  overestimates  the  length  of  the  transition  region 
and  does  not  give  a good  estimate  of  the  gradient  of  heating  in  this  region. 

Of  the  prediction  methods  considered,  the  method  of  Chen  and  Thyson  was  derived 
from  more  basic  considerations  than  the  other  methods  and  rests  upon  a broader 
base  of  experimental  data.  It  is  felt  that  until  more  extensive  comparisons  are 
made  with  experimental  data,  this  should  be  the  preferred  method  for  predicting 
transitional  heating. 

3.  HEAT  CONDUCTION  CALCULATIONS 


As  stated  in  the  Introduction,  the  present  efforts  also  include  the 
development  of  new  technology  for  aerodynamic  heating  calculations.  In  this 
section  a brief  description  of  a new  integral  method  will  be  presented.  The 
method  is  currently  under  development  into  a practical  tool  for  solving  boundary 
layer  flow  and  transient  heat  conduction  problems  alike. 

The  basic  idea  of  the  method  will  first  be  introduced  by  way  of  a simple 
example  of  an  ordinary  differential  equation.  The  characteristic  features  of  the 
method  are  explicitly  demonstrated  in  the  solution  process  of  this  model  problem 
along  with  the  principal  merits  of  the  method.  A class  of  one-dimensional 
transient  ablation  problems  is  then  solved  by  the  present  method  as  an  example 
of  application.  Many  details  of  the  method  are  omitted  from  this  paper  to 
conserve  space,  but  appear  in  Refs.  (20,  21,  22). 

3.1  Basic  Idea  of  the  Integral  Method — A Model  Example 

The  basic  idea  of  the  new  integral  method  lies  in  the  use  of  the  integrated 
version  of  the  governing  differential  equation  as  an  expression  for  the  boundary 
derivatives,  after  an  approximate  (guessed)  solution  is  substituted  for  the 
unknown  exact  solution  in  the  integrands.  In  physical  applications,  the  boundary 
derivatives  are  often  related  to  the  important  quantities  of  surface  flux,  e.g., 
skin  friction  (momentum  flux),  heating  rates  (heat  flux),  etc.  on  an  aerodynamic 
vehicle.  Their  accurate  and  efficient  predictions  are  of  critical  importance  to 
the  design  and  performance  analysis  of  the  vehicle. 


20.  Zien,  T.  F.,  "Approximate  Calculation  of  Transient  Heat  Conduction,"  AIAA  J . , 
Vol.  14,  No.  3,  Mar  1976,  pp.  404-406. 

21.  Zien,  T.  F.,  "Integral  Solutions  of  Ablation  Problems  with  Time-Dependent 
Heat  Flux,"  AIAA  Paper  78-864,  2nd  AIAA/ASME  Thermophysics  and  Heat 
Transfer  Conference,  Palo  Alto,  Calif.,  24-26  May  1978.  Also  AIAA  Journal, 
Vol.  16,  No.  12,  pp.  1287-1295,  Dec  1978. 

22.  Zien,  T.  F.,  "A  Simple  Prediction  Method  for  Viscous  Drag  and  Heating  Rates," 
Paper  presented  at  the  1978  Science  and  Engineering  Symposium  (sponsored  by 
USN/USAF) , 17-19  Oct  1978,  to  appear  in  the  Symposium  Proceedings. 
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The  basic  idea  will  here  be  illustrated  in  terms  of  a simple  boundary  value 
problem  of  ordinary  differential  equation.  While  the  model  problem  deos  not 
simulate  exactly  the  mathematical  structure  of  the  physical  problems  intended 
for  the  application  of  the  method,  it  does  serve  the  purpose  of  exhibiting  the 
general  ideas  in  an  elementary  fashion  for  easy  understanding. 


Let  us  consider  the  following  boundary  value  problem  for  the  ordinary 
differential  equation: 


0 < x < 1 


dx 

y(0)  = 0 

y(D  = 1 

The  exact  solution  is  easily  obtained  as 

sinhx 
y ~ sinhl 

whereby  the  exact  boundary  derivatives  are  calculated  as 

- 0.8059;  (^1,  = 1.3130 


(t)o  - °-8059'  (S)l  ’ !•' 


(3.1) 

(3.2a) 

(3.2b) 

(3.3a) 

(3.3b) 


These  exact  solutions  will  be  used  as  a standard  to  determine  the  accuracy 
of  the  approximate  integral  solutions  to  be  obtained. 


One  way  to  construct  approximate  solutions  to  the  system,  Eqs.  (3.1)  and 
(3.2),  is  to  use  polynomials  which  satisfy  the  boundary  conditions,  Eqs.  (3.2a,b). 
Two  such  possible  solutions  are  given  below, 

f^  = Ax  + (1  - A)x2  (3.4a) 

f2  = Ax2  + (1  - A)x3  (3.4b) 

where  A is  a constant  referred  to  as  the  profile  parameter.  It  is  to  be 
determined  by  requiring  the  approximate  solutions  to  satisfy  the  integrated 
differential  equation. 


The  integral  of  Eq.  (3.1)  has  the  form 

(£)i  - (S)o  - / 

which  will  be  used  in  the  present  method  as  the  expression  of  che  boundary 
derivatives  when  f is  substituted  for  y in  the  integrand.  In  the  elementary 
version  of  the  present  method,  we  assume  that  one  boundary  derivative,  say, 
(dy/dx)^  is  determined  by  simply  using 


(3.6) 


•**•«%>  V**SS«^ 
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Eq.  (3.5)  is  then  used  in  conjunction  with  an  auxiliary  equation  for  the 
determination  of  the  two  unknowns,  (dy/dx)Q  and  A.  The  auxiliary  equation  may 
be  generated  from  a x-moment  integral  of  Eq.  (3.1)  i.e., 

.1  .2..  .1 


/ 

0 


xpdx  = 
dx 


xydx 


(3.7) 


0 


Substitution  of  f for  y in  the  integrands  of  Eq.  (3.5)  and  Eq.  (3.7)  (using 
the  boundary  conditions,  Eqs.  (3.2),  and  assuming  Eq.  (3.6)  leads  to,  respectively, 
the  following  equations: 


and 


Eqs.  (3.8)  and  (3.9)  then  determine  the  quantities  of  A and  (dy/dx)^.  The 
results  are  given  below: 


/ A g\  / J.A 

(dxjl  " (dxjo  = / fdx 

(3.8) 

| 

\ 

/ A f\  r 

i 

(dA  - 1 = / *fdx 

(3.9) 

\dx/1  0 

} 

(i)  f = fi:  A = 9/13,  (dy/dx)0  = 0.8589  (+0.9%) 

(dy/dx)x  = (dfj/dx)!  = 2 - A = 1.3077  (-0.4%) 

(ii)  f = f2:  A = 12/7  = 1.7143,  (dy/dx)0  = 0.8929  (+4.9%) 

(dy/dx)1  = (df2/dx)x  = 3 - A = 1.2857  (-2.1%) 

If  in  the  solution  process  we  choose  to  base  the  other  boundary  derivative 
on  the  approximate  profile,  i.e.,  (dy/dx)g  = (df/dx)g,  then  Eqs.  (3.5)  and  (3.7) 
combine  to  determine  A and  (dy/dx)^.  The  results  are,  for  f = f^, 

A * 11/13  = 0.8462,  (dy/dx)x  = 7A/6  + 1/3  = 1.3205  (+0.6%) 

(dy/dx)Q  = (df1/dx)0  - A - 0.8462  (-0.6%). 


The  percentage  errors  of  these  approximate  values  of  the  boundary  derivatives 
are  included  in  the  parentheses.  These  approximate  results  are  seen  to  be  very 
accurate,  even  for  the  ot/iously  erroneous  profile  of  f2  (zero  profile  slope  at 
x = 0).  Also,  they  exhibit  rather  mild  dependence  on  the  assumed  approximate 
solution. 

It  is  useful  to  compare  these  results  with  those  of  the  classical  integral 
method  in  order  to  demonstrate  the  inadequacy  of  the  older  method  and  the  need  for 
its  improvement.  This  will  now  be  carried  out  in  the  following. 


In  the  classical  method,  the  integrated  differential  equation,  Eq.  (3.5),  is 
the  only  equation  for  the  determination  of  the  approximate  profile,  with  both 
boundary  derivatives  taken  directly  from  the  approximate  solution,  i.e.,  (dy/dx)g 
= (df/dx)g,  (dy/dx)}  = (df/dx)j.  The  classical  integral  solutions  are  given  below: 


(i)  f = fr- 

A = 10/13  - 0.7692,  (dy/dx)Q  = 0.7692  (-9.6%) 
(dy/dx)x  = 1.2308  (-6.3%) 

(ii)  f = f2: 

A = 2.5385,  (dy/dx)0  = 0 (-100%) 

(dy/dxJ-L  = 0.4615  (-65%) 

The  above  calculations  clearly  demonstrate  the  difficulties  with  the  classical 
integral  method.  It  is  generally  inaccurate,  and  very  strongly  depends  on  the 
(assumed)  approximate  profile.  The  latter  difficulty  makes  the  results  unreliable, 
as  there  exists  no  unique  way  for  choosing  an  approximate  profile  in  the  calcula- 
tion of  an  actual  physical  problem.  Note  that  the  basic  idea  of  the  classical 


’ry  -‘,y 


22 


NSWC  TR  79-21 


method  when  applied  to  boundary-layer  flow  calculations  leads  to  the  well-known 
Karman-Pohlhausen's  momentum  integral  technique. 

In  the  present  new  integral  method,  the  auxiliary  equation  can  actually  be 
generated  in  various  different  ways.  One  variant  of  the  method  is  to  use  the 
y-moment  integral  of  the  basic  differential  equation. 

In  the  y-moment  scheme,  the  auxiliary  equation  takes  the  following  form: 


)l  - fl  (1 

f)2d*  - f1  Ax 

(3.10) 

r o \d 

V o 

It  will  be  assumed  that  (dy/dx)^  = (df/dx)-^  in  the  following  solution  process. 

Solution  of  Eqs . (3.8)  and  (3.10)  gives  the  following  results  for  f = fj_: 
(the  other  negative  solution  for  A is  discarded  on  physical  grounds) 

A = 0.6826,  (dy/dx)Q  = 0.8703  (+2.3%) 

(dy/dx)1  = (df1/dx)1  = 1.3174  (+0.3%) 

Again,  the  results  are  found  to  be  very  accurate,  showing  little  variation 
from  the  previous  results  based  on  the  x-moment  scheme. 

As  a further  illustration  of  the  generalization  of  the  new  idea,  we  will 
exploit  the  combined  use  of  the  x-moment  and  the  y-moment  scheme  in  the 
approximate  solution  of  the  model  problem.  In  this  combined  scheme,  three 
equations  are  generated,  i.e.,  Eqs.  (3.5),  (3.7)  and  (3.10),  for  the  solution  of 
three  unknowns:  A,  (dy/dx)g  and  (dy/dx)^. 

For  f = fp  the  results  are 

A = 0.7727,  (dy/dx)Q  = 0.8523  (+0.2%),  (dy/dx^  = 1.3144  (+0.1%) 
which  are  practically  indistinguishable  from  exact  solutions. 


The  approximate  solutions,  f^  and  f2  are  plotted  in  Figs.  11a,  b,  c for 
both  present  method  and  the  classical  method.  It  is  seen  that  the  present 
solutions  in  the  entire  region  0 1 x t 1 do  not  show  any  substantial  improvements 
over  the  classical  ones,  although  the  boundary  derivatives  calculated  by  the 
present  method  are  far  more  accurate  than  those  by  the  classical  method.  It  is 
emphasized  here  again  that  in  the  present  method,  the  boundary  derivative  is  not 
taken  from  the  boundary  slope  of  the  assumed  profile. 

This  simple  model  problem  serves  to  bring  out  the  characteristic  features 
of  the  new  integral  method.  It  is  accurate  and  profile-insensitive,  compared 
to  the  classical  integral  method.  Its  primary  utility  is  the  calculation  of  the 
boundary  derivatives. 


In  the  following,  the  application  of  the  method  to  a class  of  idealized 
ablation  problem  will  be  presented  as  an  example  of  the  application  of  the 
method. 
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3.2  The  Ablation  Model 


The  model  used  here  is  a semi-infinite  solid  initially  in  a uniform 
temperature,  T^,  lower  than  the  phase-change  temperature  of  the  solid,  Tp.  An 
unsteady  heat  flux,  q0(t),  is  then  applied  at  the  boundary  until  the  boundary 
temperature  reaches  the  phase-change  temperature  of  the  solid.  This  period  is 
referred  to  as  the  preablation  period.  As  the  external  heating  continues, 
melting  commences  with  the  melting  front  progressing  into  the  solid,  and  this 
period  is  referred  to  as  the  ablation  period.  In  the  idealized  model,  it  is 
assumed  that  the  molten  solid  is  instantaneously  and  completely  removed  upon  its 
formation  say,  by  the  action  of  some  aerodynamic  forces,  so  that  the  melting  line 
acts  like  a new  (moving)  boundary  upon  which  the  external  heat  flux  qD(t),  acts. 
This  assumption  is  particularly  appropriate  for  the  ablation  of  subliming 
materials  such  as  comphor,  graphite,  etc.  Also,  to  simplify  the  calculations, 
the  thermophysical  properties  of  the  solid  are  assumed  constant.  This  model  is 
the  same  as  the  one  used  earlier  by  Landau23s  except  that  he  considered  only  the 
special  case  of  a constant  qQ,  and  obtained  solutions  by  entirely  numerical 
means.  The  model  is  sketched  in  Fig.  12. 

In  terms  of  this  idealized  model,  the  governing  equations  and  boundary 
conditions  are  as  follows: 


(1)  Preablation  period. 


t > t > 0, 
P 


00  > x > 0 


(3.11) 


T(x,0)  = T(°°,t)  = T^ 


(3.12a) 


- k 


(3.12b) 


where  u and  k are,  respectively,  the  (constant)  therm-sl  diffusivity  and  heat 
conductivity  of  the  solid.  Also,  tp  signifies  the  time  at  which  the  boundary 
temperature  reaches  Tp,  i.e.,  T(0,  tp)  = Tp. 

(2)  Ablation  period. 


3T 

at 


00  > t > tp,  «>  > X > X(t) 

(3.13) 

T(x,  tp)  = T(x,  t^) 

(3.14a) 

T(X,t)  = T 

P 

(3.14b) 

T(«,t)  = T 

00 

(3.14c) 

23.  Landau,  H.  G.,  "Heat  Conduction  in  a Melting  Solid,"  Quarterly  of  Applied 
Mathematics,  Vol.  8,  1950,  pp,  81-94 ♦ 
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- k (l)x-X  + «L  f ' %(t>  °-14, 

Eq.  (3.14a)  ensures  the  continuity  of  the  temperature  distribution  within 
the  solid  at  the  onset  of  ablation,  t = tp,  and  Eq.  (3.14d)  states  the  energy 
balance  across  the  ablating  front,  x = X(t).  Note  that  the  boundary  condition, 
Eq.  (3.14d),  which  relates  the  ablation  speed,  dX/dt,  to  the  temperature 
gradient,  (3T/3x)  _v,  is  the  basic  source  of  the  nonlinearity  of  the  problem. 

X”  A 

3.3  Power-Law  Boundary  Heat  Flux  - Present  Method  of  Solution 

20  21 

The  9-moment  scheme  * of  the  present  integral  method  will  be  used  in  the 
calculations  throughout  this  paper.  Briefly  speaking,  the  procedure  makes  a 
combined  use  of  the  heat  balance  integral  (the  integrated  form  of  the  heat 
equation)  and  the  integral  of  the  original  heat  equation  after  it  is  multiplied 
by  0(=  T - Tro) . A certain  approximate  temperature  profile,  f,  is  then  sub- 
stituted for  the  temperature  in  these  integrated  version  of  the  heat  equation. 
The  heat  balance  integral  based  on  the  approximate  temperature  profile  is  used 
as  the  expression  for  the  boundary  heat  flux. 

3.3.1  Preablation  Solution 


Consider  the  idealized  ablation  problem  with  qQ  = Atra,  where  A and  m are 
constants.  For  the  preablation  period,  an  exponential  profile  for  the  temperature 
excess,  0 = T - T , is  used  in  the  calculation,  i.e., 

CO 


q <5 

f = $ exp 


(3.15) 


The  profile  contains  two  parameters,  6 and  0,  and  satisfies  only  the  boundary 
of  0^  = 0.  Recall  that  the  boundary  flux  is  not  to  be  obtained  from  (3f/3x)0  in 
the  present  method20,21.  The  preablation  solutions  have  already  been  presented 
in  Ref.  (20),  and  will  be  briefly  summarized  here.  The  boundary  temperature,  Tq, 
in  dimensionless  form  is  given  as 


k(T  - T ) 

e = 2 — =Vm+5/ 4 / (m+1)  (3.16) 

0 q /at 

The  parameters  6 and  8 are 

6/v'ott  = 2/v/4m+5  (3.17a) 

0 = (m+5/4)/ (m+1)  (3.17b) 

As  is  shown  in  Ref.  (20),  the  boundary  temperature  as  given  by  Eq.  (3.16) 
agrees  better  than  1%  with  the  exact  solution  for  all  m in  the  range  0 5 m < 00 . 

Ablation  starts  when  T0  = T . From  Eq.  (3.16),  the  corresponding  time,  t , 
is  determined  as  ^ ** 
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t 

P 


- 1 

k(m  + 1) (T  - T ) m + 1/2 

p 00 

A/(m  + 5/4)ct 


(3.18a) 


The  "penetration  depth",  6,  at  t follows  from  Eq.  (3.17a), 

k(m  + 1) (T  - T )am 
6 - ■— P. 

^ A/m  + 5/4 

The  quantities  tp  and  <5p  are  conveniently  used  as  the  scales  for  time  and 
length,  respectively,  in  the  formulation  of  the  ablation  problem.  However, 
since  tp  and  6p  given  above  are  only  approximate  solutions  dependent  on  the 
method  of  solution,  it  appears  desirable  to  use  the  characteristic  time  and 
length  of  the  problem,  tc  and  &c,  as  the  scales  for  easy  determination  of  the 
absolute  accuracy  of  the  solutions.  From  a simple  dimensional  consideration, 
it  is  easily  found  that  tc  and  £c  of  the  problem  can  be  defined  by 


1 

2m  + 1 

/ ^m  + 5/4.  (3.18b) 


| 


and 


t - [k(AT)/Av^a]1/(m+1^ 
c 


(3.19a) 


£ = £T-  [k(AT)am/AJl/(2m+1) 
c c 


(3.19b) 


where  AT  = T - T . 

P 00 

Thus,  we  introduce  the  following  two  sets  of  dimensionless  variables  in  the 
ensuing  calculations: 


r = — 
* “ 6 » 


_ t 

T = — , 
P 


A - 5 
A = <$  ’ 


X = ^ 
A - 6 


(3.20a) 


and 
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(Ttp)E  = [T(m  + f)/r(^)]1/(m  + V 


(3.22) 


where  T is  the  Gamma  function.  It  can  be  easily  shown  that  t.  of  the  present 
solution  approaches  the  exact  limit  of  = 1 as  nr**'. 

The  comparison  between  the  approximate  solution  and  the  exact  solution  for 

t is  shown  in  Fig.  13.  The  present  solution  is  practically  indistinguishable 

from  the  exact  solution  in  the  entire  range  of  m,  ^ > in  > 0 in  the  figure,  the 

maximum  error  being  about  1.8%  at  m - 0.  It  is  also  interesting  to  note  the 

existence  of  a maximum  t.  near  m = 1.5. 

Ip 

3.3.2  Ablation  Solution 

The  ablation  problem  is  then  formulated  in  dimensionless  form  by  using 
dimensionless  variables  5,t,A,A  and  0.  The  normalized  temperature, 

9 =*  (T  - T ) / (T  - T ) is  also  used. 

CO  p 00 


An  exponential  profile  for  0 is  assumed. 


f = exp 


x - X(t) 
6<t) 


(3.23) 


where  X(t)  is  the  (unknown)  ablation  line  location,  and  X(tp)  = 0.  Note  that 
this  choice  of  the  temperature  ensures  the  continuity  of  the  temperature  field 
at  t ~ if  6(t)  is  assumed  to  be  continuous  at  tp. 

With  f substituting  for  0,  the  following  equations  are  easily  derived  from 
the  ablation  line  condition,  Eq.  (3.14d),  the  integration  of  Eq.  (3.13)  from 
? = A to  5 = «,  and  the  integration  of  the  8-raoment  of  Eq.  (3.13),  respectively 


- (m  -f  5/4)  -||jx  * (m  + l)rm  - v 


dA  , dA 
dT  dx 


= - (m  + 5/4)  |||. 


1 dA  , dA  , , ?//N  L 30i  . l] 

I3T  + dT-  - <”  + 5'4>  [2  a 7,h  + aj 


(3.24a) 


(3.24b) 


(3.24c) 


where  v = Qi/Cp(Tp  - T^) . Ql  is  the  latent  heat  of  ablation  per  unit  mass  and 
Cp  is  the  specific  beat  of  the  solid. 

In  the  present  procedure,  Eqs.  (3.24)  form  the  system  for  the  three  unknowns, 

A,  A and  the  heat  flux  at  the  ablation  front,  - (30/8£)^ . 

Closed-form  solution  can  be  obtained  for  the  special  case  of  m=0,  i.e. , 
the  case  of  a constant  heat  flux.  For  other  values  of  m,  solutions  are  obtained 
easily  by  numerically  integrating  an  ordinary  differential  equation  (see  Ref.  (21)). 

Typical  results  for  (m,v)  = (0,  /ir/2)  are  shown  in  Fig.  14,  compared  with 
Landau’s^  numerical  solution  and  the  classical  heat  balance  integral  (HBI) 
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solution  based  on  similar  temperature  profiles.  Note  that  the  ablation  thickness 
in  Fig.  14  is  normalized  by  its  value  at  t = «,  i.e., 


where 


X 

n 


X/X 


00 


X = T^/d+v). 

00 

It  is  clear  that  the  present  solution  is  very  accurate  and  shows  considerable 
improvements  over  the  HBI  solution. 

Admittedly,  the  example  treated  here  is  highly  idealized.  Nevertheless,  the 
essential  features  of  the  mathematical  system  of  ablation  problems  are  included 
in  the  model.  The  success  of  the  method  in  this  application  is  therefore  indica- 
tive of  its  potential  for  providing  approximate  solutions  to  realistic  conduction 
problems  in  general  and  ablation  problems  in  particular. 

4 . CONCLUDING  REMARKS 


In  the  present  paper,  an  improved  procedure  for  predicting  laminar  heating 
on  the  hemisphere-cylinder  configuration  is  proposed.  While  the  method  is  still 
empirical,  it  has  been  demonstrated  to  offer  reliable  solutions  over  a wide 
range  of  operating  conditions  and  Mach  numbers  of  interest  to  the  Navy’s  tactical 
missiles.  Inasmuch  as  the  flight  test  data  suggest  laminar  flow  over  most  part 
of  the  missile  dome,  the  improvement  should  be  considered  significant.  For 
turbulent  heating  calculations,  the  current  NWC  scheme  appears  to  be  in  reasonable 
agreement  with  the  finite-difference  calculations.  However,  since  the  exact 
solution  to  turbulent  flow  is  not  known  at  present,  the  agreement  should  not  be 
viewed  as  a validation  of  the  NWC  scheme.  The  true  confirmation  could  only  come 
from  a comparison  with  experimental  data.  In  the  transitional  region,  the  method 
based  on  the  idea  of  "spot  formation"  appears  encouraging  and  warrants  further  study. 

The  transient  conduction  calculation  presented  in  this  paper  is  an  example 
of  the  new  technology  being  developed  for  future  use  in  the  aerodynamic  heating 
studies.  As  the  same  basic  ideas  have  been  used  in  boundary  layer  calculations 
as  well,  the  encouraging  results  presented  here  offer  the  hope  of  replacing  the 
empirical  schemes  by  this  more  rational  approach  in  the  future.  A self-consistent 
procedure  for  aerodynamic  heating  calculations  will  then  become  available. 

In  the  future,  a similar  study  of  the  laminar  heating  will  be  conducted 
for  the  other  principal  dome  configurations,  such  as  the  spherical  ogive.  The 
ultimate  goal  is  to  replace  the  empirical  formulas  currently  in  use  by  some 
more  rational,  yet  still  simple,  predictive  schemes. 

The  capabilities  of  calculating  convective  heating  for  three-dimensional 
configurations  will  also  be  developed  in  the  course  of  the  study. 
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